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Linear stability analysis of an elastically anchored wing in a uniform flow is investi- 
gated both analytically and numerically. The analytical formulation explicitly takes into 
account the effect of the wake on the wing by means of Theodorsen's theory. Three dif- 
ferent parameters non-trivially rule the observed dynamics: mass density ratio between 
wing and fluid, spring elastic constant and distance between the wing center of mass and 
the spring anchor point on the wing. We found relationships between these parameters 
which rule the transition between stable equilibrium and fluttering. The shape of the 
resulting marginal curve has been successfully verified by high Reynolds number direct 
numerical simulations. Our findings are of interest in applications related to energy har- 
vesting by fluid-structure interaction, a problem which has recently attracted a great deal 
of attention. The main aim in that context is to identify the optimal physical/geometrical 
system configuration leading to large sustained motion, which is the source of energy we 
aim to extract. 
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1. Introduction 

The study of the mutual interaction between fluids and elastic o bjects is a problem 
of pa ramount importance in many fields of science and technology (jPowell &: Kenneth 
20011) . In bio- fluid mechanics, with the advent of supercomputers, it becomes a cor- 
nerstone for the quantitative understa nding of a variety o f problems ranging from blood 
pressure interaction with arterial walls (Aulisa et al. 20061 ) and the blood interaction with 
mechanical heart valv es ( de Tullio et g/lioOOf ). to animal loconiotion and self-propulsion 
(jFish fc Laudeiil2006l ) and aerodynamics of insect flights (jSand 120031 ). 
It is also a topic of growing interest in relation to the possibility of manipulating the fluid 



flow to enhance aerodynamics performances of immersed bodies ( Favier et al. 2009f ) 



Fluid-structure interaction is also a crucial aspect in the design of many engineering 
systems, e.g., aircrafts and bridges. The failure in considering the effect of oscillatory 
interactions can be catastrophic. The ultimate goal is thus to reduce at a minimum all 
sources of potentially resonant couplings between fluid and structure. 
If on the one hand the reduction of the latter mechanisms is thus a crucial need for the 
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correct project of many engineering systems, on the other hand situations exist where 
the enhancement of the same coupling mechanisms between the dynamics of flexible 
solids and surrounding air/ water flows are strongly desired and sought in order to gen- 
erate self-sustained, possibly large-amplitude, motion of the solid body. This is a typical 
requirement for energy-harvesting devices through which solid body vibration can be 
successfully converted into electrical energy. 

Among the many strategies to achieve the goal of energy extraction from solid vibra- 
tions, the idea to harvest energy from the w ind kinetic energy by flapping wings turned 
out to be particularly attractive and fruitful ( McKinnev fc DeLaurienllQSl ). It is beyond 
the scope of the present p aper to provide a review o f developments which followed from 
the seminal paper by McKinnev fc DeLauriei ( 198l[ ). There remains however much work 
ahead of us, both on the side of material science (to optimize material performances in 
term of elasticity and electrical efficiency) and on those of fluid mechanics and aeroelas- 
ticity. One of the main limitations of many of the available "pitch-and-plunge" devices 
( Hodges fc Pierc32002f) is that one or more of the active degrees of freedom of the system 
are explicitly driven (e.g., mechanically, by a motor) while the other available degrees of 
freedom are left free to interact with the flow for the final aim of producing unceasing 
oscillations. The external motor clearly reduces the net gain of the device. 



Our aim here is to propose and study from the fluid mechanic point of view a simple 
system composed by a rigid thin wing (actually, a rigid airfoil) anchored to an elastic 
spring and free to move under the action of a constant wind. By means of normal mode 
linear stability analysis and direct numerical simulations (DNS) we aim at investigating 
whether self-sustained oscillations can be obtained for suitable choices of the available 
parameters (both physical and geometrical) and in the absence of any external motor. 

The paper is organized as follow. In Sec. [2] we introduce our flapping system and the 
associated equations of motion; in Sec.[3]the linear stability analysis is presented together 
with the obtained results; in Sec. |4]we describe the numerical method we exploited both 
to verify the linear analysis predictions and to extend the study to fully nonlinear regimes. 
The last section is reserved for conclusions. 



2. The flapping system 

2.1. Equations of motion 

The system we consider consists (see Fig. [T]for a cross-sectional view) of a rigid rectan- 
gular airfoil, of lenght (chord) L, thickness 5 <Si L, and span S ^ L. An elastic spring 
connects an external flxed point to an anchor point E belonging to the longitudinal axis 
of the airfoil. Such a structure is exposed to a uniform wind U, blowing along the X axis 
(from left to right). 

We will study the cross-sectional two dimensional motion, which approximately describes 
the fully three-dimensional motion of fluid and structure for an airfoil of very high aspect 
ratio S/L. We are interested in investigating the stability of the (trivial) equilibrium con- 
figuration corresponding to the airfoil aligned along the unperturbed wind flcld U. Let 
us fix the origin of the X axis so that the coordinate of the airfoil anchor point E is in 
the origin at the equilibrium configuration {Xe = 0). Indicating by C the center of mass 
position, we define the ratio r = CE/ L, which will turn out an important parameter of 
the airfoil motion. Moreover, we will only consider situations corresponding to E on the 
left (upwind) of C. A divergence instability is indeed trivially expected in the opposite 
case, i.e. when E is on the right (downwind) of C. In this work we assume that C is in 
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Figure 1. Cross-sectional sketch of the flapping device. C denotes the wing center of mass, O 
is the wing leading edge and E is the anchor point of the elastic spring (supposed to have a 
zero rest length). The oriented curved arrows fix the convention on the positiveness of angles 
and moments. The angle a is magnified for the sake of clarity, the analysis being mostly in the 
linear regime where the wing is almost aligned along the unperturbed fiow. 



the wing midpoint (the whole analysis can be easily generalized to an arbitrary position 
of C), thus ^ r ^ 1/2. 

Let us now perturbate the equilibrium configuration and evaluate the elastic force F'^ 
and its torque. If we conveniently choose E as the pole for moments, then the elastic 
force does not produce any torque. Up to linear contributions in Taylor expansion, the 
expression for F*^ reads: 

F^ = -kxc (2.1) 



y - 'k{yc-rLe) , (2.2) 

where lowercase coordinates indicate perturbations with respect to the equilibrium con- 
figuration and k is the spring constant. 

The resulting equations of motion for the perturbations are: 



mxc = —kxc 

myc = -k{yc - r L9) ~ F^ 
1^6 = -{rLfke + rLkyc 



rLF^ 



(2.3) 
(2.4) 
(2.5) 



where Ic is the moment of inertia (around the center of mass axis), m is the wing mass, 
and Mj^ are the lift force and the lift momentum. 



2.2. Aerodynamics forces and moments 
The expression for F'" and can be obtained by treating t he point E as one does 



for the flexural point in the classical 'pitc h-and- plunge' system dHodges k, Piercell2002 ) 
and then exploiting Theodorsen's theory ( Bisplinghoff et a/]ll957l ). Because of the fact 
that the latter assumes a sinusoidal response of the system (with no damping and no 
amplification), it is fully justified to determine marginal curves in the parameter space 
separating stable regimes from unstable ones (i.e. the so-called fiutter condition). Deter- 
mining such curves is one of the main aim of the present paper. 

According to Theodorsen's theory, the following expressions arise for F'" and Mj^ (|Bispl inghoff et 
1953): 



pirb 
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where h = L/2, and and are the distances of wing points at L/4, 

L/2 and 3L/4 from the leading edge. All terms involving b are the so-called added- mass 
terms. 

The function C{q) is the Theodorsen's function t hrough which the wake effect on the 
wing motion can be explicitly taken into account (jBisplinghoff et a/.l 119571 ) . It depends 
on the reduced frequency q = ujL/(2C/), where uj is the pulsation associated to the 
system harmonic response. T he form of C{q) can be given in terms of Hankel functions 
(Cradshtevn fc Rvzhikiri980l) as: 



C{q) = 



H?{q) 



(2). 



(2.8) 



Eqs. (I2.3p - (j2.5p with the expressions for the aerodynamics force and moment (|2.6p and 
(|2.7|) and the relationships xc = xe, yc = Ue + r L9 constitute our equations of motion 
for Xc, yc and 6 under the constraint of harmonic behavior, as dictated by Theodorsen's 
theory. 

Some comments on the structure of the equations above are worth discussing. The equa- 
tion for Xc is decoupled from the others; this variable can be thus ig nored. With respect 
to the classical 'pitch-and-plunge' problem (jHodges fc Piercd 120021 ) here yc and 6 are 
two-way coupled in a stronger way, a consequence of the elastic force acting on the wing 
which jointly depends on both yc and 9. This is expected to cause interesting and peculiar 
behaviors. 



3. Linear stability analysis 

Let us now investigate the normal mode linear stability analysis of the system. This 
can be easily done by normal mode decomposition: 

yc{t) = Ye"^' + C.C e{t) = Q e"^' + c.c. , 

where c.c. stands for complex conjugate. 

Plugging the expressions above into Eqs. ()2.4p and ()2.5|) we obtain an algebraic linear 
system for the two (complex) amplitudes Y and 0. The condition of having nontrivial 
solutions provides a set of two real equations for the real and immaginary parts of uj. 
Focusing on the determination of marginal curves associated to flutter conditions, and 
thus assuming the imaginary part of w to be zero, one of the two equations furnishes the 
relationship between parameters identifying the marginal curve, e.g, p as a function of k 
and r. 

Due to the presence of Theodorsen's function the determination of the marginal curve 
requires the exploitation of standard searching methods to find zeros of the associated 
(non polynomial) set of equations. This is however a task which can be easily achieved 
by means of a classical Raphson-Newton method. 

As customary, it is convenient to pass to a dimensionless form of all variables. For this 
purpose, let us define the following dimensionless quantities: 

j_U yc k Pw S 

t-^t —Tri^^ ^^Pw , 3.1 

L L ptV^ pL 
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Figure 2. The marginal curve for r = 0.5 (left, corresponding to the leading edge anchor point) 
and r = 0.26 (right). Circles refer to points analyzed through DNS. For filled circles DNS predict 
unstable motion; for open circles they predict stability. For comparison, dotted lines represent 
the marginal curve under the quasi-steady hypothesis. To emphasize variations along the vertical 
axis, a log scale is used for p„ . 

SO that the reduced frequency is simply q = uj /2. To define the dimensionless parameters 
above, we have defined by p the fluid density and written the wing mass to as to = p-uj5 L, 
pw being the wing density. 

The behavior of the marginal curve for r — 0.5 and r ~ 0.26 is reported in Fig. [2l 
Circles refer to points analyzed by DNS for Re=10000 and Re=60000. Their discussion 
is postponed to Sec. 14.21 Dotted lines represent the marginal curves under the quasi- 
steady hypothesis (i.e. when C{q) = 1 and the equations defining marginal curves are 
polynomials). This corresponds to neglecting the effect of the wake on the motion of the 
wing. Some points are worth discussing. The first is on the shape of the marginal curves 
obtained. They are quite structured and the relationship between the dimensionless pa- 
rameters associated to the fluttering condition is far from being trivial. The pronounced 
'nose' for r = 0.5 around pw ~ 10, 1/fc ~ 0.2 is an example. This peculiar shape leads 
to interesting consequence on the fluttering condition in terms, e.g., of the wing mass: 
there are indeed three different values for the wing mass (actually for the dimensionless 
quantity S / {p L)) leading to fluttering for a fixed spring constant (actually for a fixed 
dimensionless parameter k/{pU'^)). It turns out that, for all spring constants and wind 
speed, it is impossible to sustain fiapping if pyj6/{pL) < 0.59. The latter threshold has 
been determined by the quasi-steady assumption that tends to coincide (see dotted line 
in Fig. [2]) with Theodorsen's marginal curve for small spring constants. On the other 
hand, for very large value of p^ S/{pL) (order of 300 for r = 0.5; a much larger value for 
r = 0.26) fluttering is expected at all velocities. 

Interestingly, the conclusions above are very sensitive to the value of r, i.e., to the relative 
distance between the wing center of mass and the anchor point (see right frame of Fig. [21). 

To better understand the sensitivity to r we report in Fig. [3] the marginal curves in 
the plane pu,-r for two fixed values of k. The non-trivial role played by the couple p^-r 
in the emergence of stable/unstable regimes can be clearly detected from the figure. 

4. Direct numerical simulations 

The aim of this section is to give a brief description of the numerical method used 
to numerically approximate the laminar incompressible Navier-Stokes equations (so in 
essence we are doing direct numerical simulations or DNS) and verify the predictions 
obtained from the normal mode linear stability analysis with the solutions obtained from 
the DNS simulations. 



6 



A. Orchini, A. Mazzino, J. Guerrero, R. Festa, C. Boragno 




10"' lo" lo' 10- 10^ 10"' lo" lo' 10^ lo' 

Pw Pw 



Figure 3. The marginal curve for k — 1/0.2 (left) and k — 1/0.17 (right) in the pm-r space. 
'Ui' and 'f/2' denote Unstable regions with one and two unstable normal modes, respectively; 
'S' stands for Stable. For filled circles DNS predict unstable motion; for open circles they predict 
stability. 



It is beyond the scope of the present paper to fully analyze the shape and position of 
the obtained marginal curves in the whole parameter space k-pw This issue is indeed 
quite expensive from the numerical viewpoint and it will be eventually postponed to a 
possible future analysis. Rather, our aim is to corroborate our former, and to some extent 
surprising, conclusions on the nontrivial shape of the marginal curves. 

4.1. The numerical method 

Hereupon, we briefly outline the solution methodology used to solve the governing equa- 
tions on moving overlapping structured grids. The complete descri ption of th e nume rical 
method and grid ding r netho dology can be found in the papers by iHenshawl (|1994l ) and 
Ichesshire fc Henshawl 

In this manuscript, we numerically approximate the laminar incompressible Navier- 
Stokes equat ions by using the velo city-pressure formulati on or pr essure-Poisson equation 
(PPE) (HenshawUlOoS ISani et a/][2006.; Petersson.200]]: lGresho|[l99h ). 



— + udu 

ot 

— 4 
P 



h I'd u 



for x€V, t> 0, 



daU ■ dua = for X G 2?, ^ > 0, 



(4.1) 
(4.2) 



with the following boundary and initial conditions 

B{u,p) =0 for 

for 



d u =0 
u(x,0) = uo (x) 



for 



xedV, t > 0, (4.3) 
X e 02?, t > 0, (4.4) 
xeV. (4.5) 

— {x, y) (for N 2 where N is the 



In this initial-boundary- value problem (IB VP), x 
number of space dimensions) is the physical space coordinates, I? is a bounded domain 
in K^, dV is the boundary of the domain T), t is the physical time, u = (u, v) is the veloc- 
ity field, p is the pressure, v is the kinematic viscosity, p is the density, i? is a boundary 
operator initial data. 

and Uq are the initial conditions. Hence, we look for an approximate numerical solution 
of equations (|4.ip and (j4.2p in a given domain P, with prescribed boundary conditions 
and given initial conditions (equations (I4.3l) - (|4.5p ). Equations (j4.ip - (|4.5p are solved in 
logical ly rectangular grids in the t ra nsformed comput a tional s pace C = C(^,ri , t) (re - 
fer to Chesshire &: Henshawl (|l990l ): lOrikakis fc Rided (|2004[ l: iGuerrerol (|2009[ l2010f l: 
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[Henshawl ( 19941 ): Vinokuij ( 1974 ) for a detailed derivation), using second-order centred 



finite-difference approximations on structured overlapping grids. 

In general, the motion of the component grids Qg of an overlapping grid system G = {Gg}, 
may be an user-defined time dependent function, may obey the Newton-Euler equations 
for the case of rigid body motion or may correspond to the boundary nodes displace- 
ment in response to the forces exerted by the fluid pressure for the case of fluid-structure 
interaction problems. For moving overlapping grids, Eqs. (|4.ip - (l4.2p are expressed in a 
reference frame moving with the component grid as follows, 



du 



u G] d 

d^p 



u = -— + i'd^u for xeD, t>0, (4.6) 
P 



daU -dua^O for x e 2?, t> 0, (4.7) 



where G is the rate of change of the position of a given set of grid points x^ of a component 
grid Qg in the physical space (grid velocity). It is important to mention that the new 
governing equations expressed in the moving reference frame must be accompanied by 
proper boundary conditions. For a moving body with a corresponding moving no-slip 
wall, only one constraint may be applied and this corresponds to the velocity on the 
wall, such as 

U {x^\wall,t) = G (-K^lwalht) , whcrC X^lwall & dV^all {t) ■ (4.8) 

Finally, in order to keep the solution of the pressure equation decoupled from the solution 
of the velocity components, we choose a time stepping scheme for the velocity co mponents 
that only involves the pressure fro m the previous time steps (split-step scheme) ()llenshaw 



19941: iHenshaw fc Peterss"o3l2003t ) . Then, the discretized equations are integrated in time 



using a semi-implicit multi-step method, that uses a Crank-Nicolson scheme for the 
viscous terms and a second-order Adams-Bashforth/ Adams- Moulton predictor-corrector 
approach for the convective terms and pressure. This solution method yields a stable 
second-order accurate in space and time numerical scheme on moving overlapping struc- 
tured grids. 

To assemble the overlapping grid system G and solve the laminar incompressible Navier- 
Stokes equations in their velocity-pressure formulation, the Overture^ framework is used. 
The large sparse non-linear system of equations arising from the discretization of the lam- 
inar incompressible Navier-Stokes equations is solved using the PETScU library, which 
was interfaced with Overture. The system of non-linear equations is then solved using a 
Newton-Krylov iterative method, in combination with a suitable preconditioner. 

As a final remark, for the range of Reynolds numbers considered in this study, no 
turbulence model or sub-grid scale model is used for the numerical simulations conducted, 
so in essence two-dimensional direct-numerical simulations (DNS) are performed. Despite 
the fact that the grid does not resolve the smallest scales in the far-wake region (where 
small scale activity is however strongly reduced) , the grid resolution is always more than 
adequate. It is so also in the near- wake flow field where the resolution is sufficiently high 
to properly resolve the small scales so that accurate instantaneous drag and lift forces 
are computed. 



f https:/ /computation. llnl.gov/casc/Overture/ 
X http:/ /www-unix. mcs.anl.gov/petsc/petsc-as/ 
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Figure 4. Time behavior of the pitching angle 9 (deg), on the left, and of the j/-coordinate of 
the wing leading edge, on the right, for the unstable case with 1/k = 0.35 and = 3.5 and 
r — 1/2. All variables are made dimensionless by p.ip . 
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Figure 5. The same as in Fig. |4]but for the stable case 1/k — 0.1 and pu 



20. 



4.2. Numerical settings and results 

We confine our attention to the marginal curve of Fig. [2] (left), performing direct numer- 
ical simulations for different values of {p^, 1/^), for Reynolds numbers between 10000 
and 60000. We almost never place ourselves too close to the predicted marginal line to 
avoid a very slow increase/decrease of the initial perturbation and thus to avoid long and 
expensive simulations. As far as initial perturbations are concerned, we firstly consider 
the one dimensional problem (along x) and we reach the equilibrium position under that 
condition. At this stage we impose a perturbation on the angle 9, of the order of 10"^ deg, 
and leave the system free to evolve. 

Filled circles in Fig. [2] (left) are associated to unstable behavior; stability is represented 
by open circles. The agreement with our predictions based on Theodorsen's theory is 
good. 

Typical time behavior for the angle 9 and the leading edge vertical coordinate yo are 
reported in Fig. H] for — 3.5 and 1/k — 0.35 (unstable case, see also the movie) and 
in Fig. [5] for = 20 and 1/k = 0.1 (stable case, see also the movie). Notice how the 
unstable case reaches a fully nonlinear stage characterized by a sustained flapping with 
a y excursion of the wing leading edge of order one (with lengths normalized by the 
wing chord) . This is a very promising result in view of applications related to the energy 
harvesting where elastomeric capacitors need to be stretched/compressed much, in order 
to produce reasonable amounts of electric energy. 

The example shown in Fig. |4] tells us that an infinitesimal perturbation yields the sys- 
tem to a finite- amplitude limit cycle associated to unceasing flapping. A natural question 
which arises is on whether unceasing flapping can be observed if a finite size perturbation 
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Figure 6. Time behavior of the pitching angle 9 (deg), on the left, and of the j/-coordinate of 
the wing leading edge, on the right, for the linearly stable case with 1/k — 0.1, = 20 and 
r = 1/2. All variables are made dimensionless by (|3.1|) . The initial perturbation is here of finite 
size. 

is applied to configurations which are stable with respect to small perturbations. As an 
example we address this question for the couple of parameters pw — 20 and 1/fc = 0.1. 
The initial angle is now 90 deg. The resulting time behaviors for 9 and yo are reported 
in Fig. ini We do not detect any tendency of relaxation toward the aligned configuration, 
a fingerprint of a subcritical instability. 



5. Conclusions 

A simple flapping system has been presented and investigated both analytically and 
numerically. Three dimensionless parameters come into play and their role on the result- 
ing flapping states appear to be highly non trivial. This property has been confirmed by 
numerical simulations we carried out in the range of Reynolds numbers between 10000 
and 60000. 

Both supercritical instabilities and subcritical instabilities have been found in our system. 

Our findings have direct applicability to the energy harvesting problem by fluid- 
structure interaction (jBoragno et aZ.1 120121 ). We have found sets of parameters leading 
to unceasing flapping for which the excursion of the leading edge is of the same order of 
the wing chord. Once our spring is replaced by an elastomeric capacitor, the observed 
large amplitude of the wing oscillations is a necessary condition for extracting a reason- 
able amount of energy from the device. 
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